A NUMERICAL MINIMIZATION SCHEME FOR THE COMPLEX 

HELMHOLTZ EQUATION 

RUSSELL B. RICHINS AND DAVID C. DOBSON 



C ^ ^ Abstract. We use the work of Milton, Seppecher, and Bouchittc on variational principles for 

^Nj waves in lossy media to formulate a finite element method for solving the complex Helmholtz 

I equation that is based entirely on minimization. In particular, this method results in a finite 

^ element matrix that is symmetric positive-definite and therefore simple iterative descent methods 

^~5 and preconditioning can be used to solve the resulting system of equations. We also derive an 

fy. error bound for the method and illustrate the method with numerical experiments. 
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1. Introduction 



(N 

< 

^^ Many systems that result in steady-state oscillations can be modeled with the Helmholtz equation, 

• but of particular interest are acoustic waves and transverse electric or transverse magnetic electro- 

'^j magnetic waves in inhomogeneous media. In each of these situations, the equation of interest can 

^ be expressed as 

I 1 (1) -V-p^^VP P = 

K 

for appropriate choices of the complex-valued, spatially dependent material parameters p and k, 
where w > is the frequency. The classical methods of deriving a weak form for this equation (with 
Dirichlet boundary conditions, for example) result in the variational equation 



K 



,2 



"^VP- Vu Pu 



dx = o, VueiJo^(r), 



f~^ which corresponds to a stationary principle, but not a minimization principle. In [7], Milton, 

^^ Seppecher, and Bouchitte expand upon the work of Cherkaev and Gibiansky [3] for the conductivity 

^~~^ equation to derive variational principles for (fTl) (as a special case of the more general equations of 

^ elasticity and electromagnetism) that are true minimization principles, provided the media are 

lossy. The minimization functional corresponds physically to dissipated energy in the system, and 

is valid even for arbitrarily small coefficients of loss. While the framework presented in [7] results in 



^ nonstandard boundary conditions, Milton and Willis extend the principles to handle the classical 

Dirichlet and Neumann boundary conditions in [8j. 

In this paper we apply the finite element method to these minimization principles and thereby 
develop a numerical algorithm for solving (fTl) that can take advantage of the many efficient methods 
available for solving a symmetric, positive-definite system of linear equations. The outline of the 
paper is as follows. Sections [2] and [3] review the general variational formulation and boundary 
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conditions introduced by Milton, Seppecher, Bouchitte, and Willis. In Section|4J we derive an error 
bound on certain finite element discretizations of the variational principle. In Sections [5] and l6J we 
describe a straightforward implementation of the finite element method on a square domain, with 
Dirichlet boundary conditions. In Section [71 we suggest a preconditioner for solving the resulting 
symmetric positive definite linear system via the preconditioned conjugate gradient method, and 
find conditions on the material coefficients under which we expect the best conditioning. Section [8] 
describes the results of some simple numerical experiments, and illustrates numerical convergence 
consistent with the error bounds from Section |4J Finally, in Section |9J we extend the method to 
handle Robin boundary conditions, and present some associated numerical examples. 

2. Variational Formulation 



Our model problem is 








(2) 


1 -V-p-^VP- 
1 P-f 


hi 


= inP 

onar 



where F is an open, bounded subset of M"^ (d = 2 or 3) with smooth boundary. For acoustic waves, 
p is the density, n is the bulk modulus, w is the frequency, and P is the pressure. Here p, k, and P 
are all complex. In this section, we focus on Dirichlet boundary conditions for simplicity; Neumann 
conditions can be handled similarly. In [7], it is shown in detail how this and other problems can 
be formulated as a minimization. What follows is a brief outline of the general framework. 

Let J-{x) and Q{x) be complex- valued fields of the form 



J J \ 9 

where F, G : F — > C'* and /, g : F — > C Suppose there exists a complex- valued potential u such 
that 



and that Q satisfies 



J" = nw 

u 



h + ug^Q, 



Z = I r^T 



where UQ := —V -G + g and /i is a source term. Suppose also that J- and G satisfy the constitutive 
relation 

(3) Gix) = Z{x)T{x) 

where Z has the form 

L K 
K^ M 

Then the constitutive relation along with the differential constraints imply 

(4) h + U{Z n u) = or V • (LVu + Ku) = h + K^Wu + Mu 

Let n be the unit outward normal on dV . It is shown in [7] that if we are given Uq and Gg • n and 
we specify 

(5) u' = u'q and G' • n = Gq • n on dV 
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(herein ' denotes the real part of a complex quantity and " the imaginary part), then the solution 
to Q satisfying the boundary conditions ([5]) is a minimizer of the functional 



Y{u',G') = j 



-Q' ) \ -Q' 



dx. 



where 

G" \ ( Z" + Z'{Z")-^Z' Z'{Z")-^ 

V -G' -h ) ^"^'^^-y {Z")'^Z' {Z")-^ 

provided that C is positive definite. An inspection of the constitutive relation shows that C is 
positive definite as long as Z" is. Explicitly, following [7j we see that if we let J-"' and Q' be 
arbitrary, and define Q" and F" by 

f Q" \ r( T' 

which is equivalent to 



T" ) ~^\ -G' 



then 



g' = Z'F' - Z"F" 
g" = Z'T" + Z"T' ' 



'^g')-^{^g' ]=:F'-g"-F"- g' 

7II -ci\ 



= T' ■ {Z'T" + Z"T') - T" ■ {Z'T' - Z"T" 
= F' -Z''^' + F" ■Z"F". 
Therefore, C is positive definite as long as Z" is. 



In addition to the conditions 



3. Boundary Conditions 



u' — u'q and G' ■ n = G'^- n aiv dV 



we can also solve the problem for u' and G' with the boundary conditions 
(6) u" = u'^ and G" ■ n = G'^ ■ n on dV, 

u' — u'q and u" — Uq on dV , 
or G" • n = Go • n and G" ■ n = Gq • n on dT. 

The correct variational principles for the last two sets of boundary conditions can be deduced from 
the formulations for the first two. The second boundary condition above is a condition on the dual 
(imaginary) variables u" and G", and therefore it may be enforced through boundary integrals, as 
follows. 

For simplicity, suppose h = Q. Let s G H^(T) and T e 7J(div, F). If u and G are such that the 
differential constraints and constitutive relation are satisfied, then multiplying by s and integrating, 
we get 

= / {ugy's dx= f [(-V • G" + g")s] dx= f [(-V • G" + g")s - T ■ {Vu" - Vu")] dx. 
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Integrating by parts, we find 

/ [G" ■ Vs + g"s - T ■ Vu" - u"V -T] dx ^ f [sG" ■ n - u"T ■ n] dS. 
Jr JdT 

Let T = (T, V • T)^. The left-hand side above can be re- written as 
The corresponding functional for the boundary condition 



Uq and G" ■ n = Gg ■ n on dV 



is then 



Y{u\ G')= j ( ^f Yc(^^']dx + 2f [u'^G' ■ n - u'G'^ ■ n] dS. 
To solve the PDE with the Dirichlet boundary conditions we minimize the functional 



Y{u',G')^Y{u' + u'o,G') + 2 u^'G' 

Jdr 

over u' G Hq{T) and G' G -ff (div, F). To solve the PDE with the Neumann boundary conditions we 
minimize the functional 

Y{u', G') = Y{u', G' + G'o)-2 f u'G'^ ■ n dS 

Jdr 

overw' e H^{r) and G' G Ho{diy,T) = {v G i7(div,r) : {vn,w) ^ Oy w e H^iT)} (see %). 

4. Error Bound 

In this section we give a bound on the error incurred by solving any of the minimization problems 
above over a finite dimensional subspace of the specified Sobolev spaces. We will give a more 
detailed account of exactly what the finite dimensional space looks like later on; in this section all 
that will matter is the highest degree of polynomials that the finite dimensional space contains. We 
will use the Bramble-Hilbert lemma to give a bound on the error. 

Here we will drop the primes used to denote real and imaginary parts. Note that what follows 
applies to any of the boundary value problems discussed previously, since the bounds depend only 
on the corresponding bilinear form. Throughout this section, G is a constant independent of the 
solution (P, v) and the grid spacing h. 

4.1. Bilinear Form. Define the bilinear form B by 

(7) B{P,v; s,T) = Jf^^Ycf^^) dx, 

Where, as before, T = Hu, Q — (G, V-G)"^, and S and T are generated from test function s G H^iX) 
and T G iJ(div, F) in the same fashion. Assume that there exist constants 71,72 > such that 
C > 72/ and that [C{x)]ij < 71 for a.e. x G F. Let V = Hl{V) x iJ(div, F), endowed with the 
norm 

||(z.,G)||v = (||^i||l,.(r) + l|G||^(div,r))^- 
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Then it follows immediately from (It]) that 

(8) B{u,G;s,T) < C-/i\\{u,G)\\v\\{s,T)\\v 
and 

(9) B{u,G;u,G)>j2Uu,G)\\l. 

4.2. Minimization Inequality. Define an energy by 

f{s,T)^^B{s,T;s,T)-Fis,T), 

where F : H^{T) x iJ(div, F) — > M (in practice, F is usually composed of terms resulting from 
an inhomogeneous term and enforcement of the desired boundary conditions). If (u, G) is the 
minimizcr of the energy, then this pair must satisfy the Euler-Lagrange equation 

B{u,G;s,T)^F{s,T) Vs e -ff(J(r), V T G il(div, T), 

so that 

/(s, T) = f{u, G) + ^B{u -s,G-T;u-s,G-T) V seH^{r)\fT e H{dW, T). 

Consider a finite dimensional subspace Vn — Vjvi x Vjv2 of V, where Vni is a finite dimensional 
subspace of H^{r) and Vn2 is a finite dimensional subspace of iJ(div, F). If {un,Gn) is such 
that 



then 



/(ujv,Giv)= min /(s,r), 
{s,T)eVN 



[B{u — UN, G — Gn; u — Wat, G — Gn)] ^ — niin [B{u ~ s,G — T;u — s,G ~ T)] ^ 

{s,T)eVN 



Inequalities Q and ([9| imply that 

Vl2\\{s,T)\\v < VBis,T;s,T) < C^i\\{s,T)\\v V {s,T) e V, 
so we have 



(10) Vl2\\iu,G)-{uN,GN)\\v< , min C^\\iu,G) ~ is,T)\\v. 

Let Fi be the orthogonal projection from H^{T) onto V^i. Since Fi is an orthogonal projection, 
it has \\Fi\\B(H^{r),H^{r)) = 1: where B{H^{T),H^(r)) is the set of bounded linear functions from 
H^{r) to H^{T). Also, define an operator F2 : iJ (div, F) -^ Vn2 by the solution of the variational 
inequality 

{F2G, Q - -F2G)^2(r Rd) > {G,Q - -F2G)^2(r,E<i) V Q e Eq, 

over the set Eq = {v £ Vn2 '■ ||V • v\\l2(y) < ||V • G||L2(r)}, which is a closed, convex subset of 
L'^iT^R'^). We then have 

ll-^2G||^2(r,Rd) = (-F'2G,F2G)^2(r,Ed) < (G, -F2G)^2(r,Rd) < ||G||L2(r,Md)||-F2G||i2(r,E<i). 



If we take s = Fiu and T — F2G in ( 10 ), then we have 



(11) \\{u,G)-iuN,GN)\\v<C\\iu-Fiu,G-F2G)\\v. 
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4.2.1. Seminorm bounds. We will discretize the domain F by by subdividing it into smaller regions, 
each of which can be seen as a suitable shifting and scaling of a reference element. More precisely, 
if e is our reference element, there exist afhne changes of variables Fi{x) = Bx + xi such that 
Fi{e) = ej, where e/ is the lih element (subdivision) in the finite element decomposition of T. In 
the case of rectangular elements in R'^, for example, we can take e = (0, 1)"^, and then we have 
B = hid- In this section a hat will denote the corresponding function defined over the reference 
element. 



Let 

(12) [uMs=Y. 

\a\ = s 

where for vector functions we define 



D'^u ■ D°'w dx and 



D^w 



D°'W2 



V D"wd J 



From 


m 


we 


get the 


inequalities 




(13) 










2 



2 \W\s^ei < \w\., < ch" 



'2 |u;| 



h'+i\V ■ q\s,e, <\V-q\s<h'^i\V-ql,e, 
for B — h, scalar functions w, and vector functions g, where w — w o F^^ and q = q o F^^ and 



denotes ( 12 ) with e; in place of e. 



We now recall the following lemma from [T, which will be used in what follows. 

Lemma 1 (Bramble-Hilbert Lemma). For some region fi C M^ and some integer k > —1, let there 
be given a bounded linear functional 

satisfying |/(w)| < <5||u||jLffc+i(Q) for all u E H''^^{n) for some 6 independent of u. Suppose that 
f{u) = for all u e Pfc(f2). Then there exists a constant C , dependent only on 57 such that 



\f{u)\<C5\u\k+u ueH'+'m 



V -q e H^{i)}. For fixed 



Let us suppose that P e H''+'^{e) and v e iJ^(div, e) = {q e W {e,\ 
elements w G H'^{e) and Q e 7J^(div, e) define the functional 

f^{u) = [u~F^uMs, /2(G) = [G-F2G,g]o, /3(V • G) = [V • G - V • F2G, V • Q]o, 
where s = or s = 1. Then, since 

\fi{u)\ <\u~ Fiul\w\s < {\uU + \Fiu\s)\w\s < (|lu|lHi(r) 
< '2\\u\\m(r)\w\s < 2||M||fffc+i(r)|u'U, 
|/2(G)| <\G- F2G|o|g|o < (|G|o + |F2G|o)|g|o = {\\G\\mrM- 
< 2||G||2,2(r)^Rd)|(3|o < l|G||fl-3(r,B<')l'3lo, 
|/3(V • G)| < |V • G - V • i^2G|o|V • Qlo < (|V • G|o + |V • i^2G|o)|V • Q|o 



\Fiu\\Hi{r))\w\s 
' ll^2G||i2(r^Rd))|Q|o 
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- (I|V • G|U2(r) + ||V • F2G|U2(p))|V • Q|o < 2||V • G|U2(r)|V • Q|o < 2||V • G||H.(r)|V • Q|o, 
and Fiu — u for polynomials in V^i and F2G — G for vectors of polynomials from Vi^2i we can 
apply the Bramble-Hilbert lemma to find that there exists a constant such that 

\h{u)\ < C\wU\u\k+i, |/2(G)| < C\Q\o\G\„ \MV ■ G)\ < C\V ■ Q|o|V • G|„ 

as long as k and j are small enough so that all polynomials of degree less than or equal to k are 
contained in the span of the basis functions representing ii and all polynomials of degree less than or 
equal to j are contained in the span of the basis functions representing G. By choosing w — u — Fiu 
and Q = G — F2G, we find that 

1^ - FiuU < C\u\k+i, \G - F^Glo < C\G\,, |V • G - V • F2GI0 < G|V • G|,. 



Employing (13), we find that ioi h < 1, 

\u - FiuU^,, < Ghi-'\u - F,uU < Ch-2-'\u\k+i < Ch''~'+^\u\k+i,e„ 
|G-i^2G|o,e, < h^-^G-F2G\o < h^-^G\G\, < G/i^'|GU, 

|V • G - V • i^2G|o,e, </i"^|V-G-V-F2G|o </i"^G|V-G|j < Ch^V ■ G\j,er 
Returning to inequality (11), we have 

\\{u,G) - {un,GnWv < C\\{u,G) - (^iw,F2G)||^ 

= G^ [|u - F,u\l,^ +\u- F,u\l,^ + \G~ F2G|g,,, + |V • « - V • F^GlU 

<gy: [h^'^MUe. + /^^l^il+M, + h'^G\l^ + h'^v ■ G\ij 

I 

<G{h''\u\l^,^r + h''{\G\lr + \\7-G\lr)). 

Let Pfe(r) denote all polynomials of degree less than or equal to k on F. We have now proved 

Theorem 1. If the solution (u, G) S H^~^^{r) x H^^^{div,T) and the finite element subspace used 
in the numerical method contains PkO^) x Pj{^)) ^ -fj(r), then there is a constant C such that the 
error satisfies 

\\{u,G) - {un,Gn)\\1 < C(/i'1«|^.+i,r + h'H\G\lr + |V • G\lr)), 
where h < 1 is the grid spacing. 

4.3. Regularity. In order for the error bound to be meaningful, we must have fc,j > 1 in Theo- 
rem [l] which means that at least 

u e H^{T) and G G iJ^(div,r). 

In the notation of the acoustic equation, if p~^ is positive definite, bounded, and G^, then classical 
elliptic regularity theory such as in [B] guarantees that P' S H^iT). Also since 

V = --p-^VP, 

we have that v E (7?^(r))^, and multiplying the acoustic equation through by —l/cu tells us 
that 

V • w = — P, 

K 
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SO V • w G ^^(r) as long as k is at least C^. 

It would be more satisfying (and useful in other contexts) to have a regularity theory derived from 
the weak form of the equations presented herein, and this is a current topic of inquiry for the 
authors. 

5. Euler-Lagrange Equation for the Model Problem 

For our model in the development of the numerical method, we will focus on the Dirichlet problem 
with functional 

f (.', G') = /^ ( ^^: ) . £ ( % )d. + 2 j^^ a . n." ,S. 

Suppose that v! and u" satisfy (ro| and (u',G") minimizes Y over all u' g Uq + -ffd(r) and G" G 
H(div, r). Then if we take any functions s e i?o (r) and T e iJ(div, T) and let T ={t,\I ■ T^ , we 
have that 

nu' + ts,G' + tT)= /f n,/ + tn.sy^/n«' + tn. 



n?i' + 1 n s 
r V -Q' - tT 



-Q' - tT 



dx 



+2 / (G" + tT) ■ nu" dS 
Jar 
has a minimum at t = 0. Therefore, 

(14) = 2/^(^^:)./:(^^)d. + 2/^^T.W'd5. 

This is the weak form of the equation that we want to solve for m'. In the case of the acoustic 
equation ([2|, we have 

u = P, L = -p-\ K ^0, M = uj^/k, h = 0, v = {~i/u})p-'^VP, G = -iujv, 

so we can rewrite 



Y{P',v") 



VF' 



n 



VF' 



ujP' 

-V • v" 



K, 



LUP' 

-V • v" 



dx 



(15) 



9r 



where r ~ —p ^, k — k ^, and 

■ r" + r'{r")-^r' r'{r")-^ 



7^ = 



(r")- 



/ k" + {k'Y/k" k'/k" 



The requirement that Z" be positive definite translates to the requirement that 
(16) p" >aT k" <-/3, a,/3>0. 



Making the substitutions in ( 14 ) for the acoustic equation, we find that the Euler-Lagrange equation 
becomes 



(17) = ^ 



VP' 



■ n 



Vs 
-uT 



ojP' 

-V • v" 



K. 



LdS 

-V-T 



dx 



LoT ■ nP" dS 



dr 



for any s G H^(r) and any T e iJ(div, T). 
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6. Implementation of the Finite Element Method 



Our goal is to test the efficacy of this new variational principle, using a simple, explicit finite 
element implementation. Let us assume that d — 2 and F = (0, 1)^. In order to find a numerical 
solution for P' , we introduce an A^ x A^ computational grid with equally spaced nodes {xj,yt) for 
t,j = I, 2, . . . , A^ and grid spacing h — 1/{N— 1). We also introduce the finite element spaces 

* = span ; M _ 1^ - ^■' 




where 



1 [i\x-Xj\,\y-yt\<h 
otherwise 



Xtjix^y) -- 

The bases of each of the finite element spaces are built from simple piecewise bilinear elements. 
We can re-index these elements with a single index by setting 



■0fc 



\x — Xj\ 



Hk = 



02k 



h 

\X-Xj\ 



\y-yt\ 



Xtj, where k^{t-2)iN~2)+j-l, k ^ 1, . . . ,{N - 2)^, 



yt\ 



I- 



\x~x^\ 



1- 



Xtj where k — {t — 1)N + j, k — 1, . 



y-yt\\ Xtj where k = {t - l)N + j, fc = 1, . 



,N' 



,7V2 



h J \ h 

We assume that our finite element solution has the form 



,Ar(JV-l) 



^R^YX=\ '^feV'fc 



^N{N-\) 



Here i/'.r is any function that satisfies the desired Dirichlct boundary condition for P' . Making this 
substitution into (17), we get 



E^vV'fc 

-CJ Y^ M\k ~^Y. Ik<p2k 



■ n 



Vs 
-ujT 



ujY.Sk'fPk 
-EPk^-^ik-Eik^' 



02k 



K, 



VV-o 




• 7^ 



Vs 

-ujT 







/c 



LOS 

-SI -T 



LOS 

-V-T 

dx — 



dx 



/ [wVV'/ • T + Lui^iV ■ T] dx, 



where we have used the divergence theorem on the boundary integral, ^i is any function on F 
satisfying the desired Dirichlet boundary condition for P", and s S H^iT), T E i?(div,F) are 
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arbitrary. In particular, this must hold when 

s = i}k, T = for fc = 1, . . . , (AT - 2)2 
s := 0, T = (/-ife for A: := 1, . . . , N{N - 1) 
s = 0, T = 02fe for k^ 1,...,N{N- 1). 

This gives rise to a system of equations of the form Aa = b, where A has the block form 

/ Ai Ai As 
(18) A=\ Ai A2 As 

V Aq As A3 

and the blocks have entries 



The right-hand side vector b is partitioned as 



b=\ b2 
bs 



where 



(6i)fc= -^ 



ib2)k= - 



[ 



UltpU \ jyf ^i'k 



^i>R 



dx 



■ n 







(20) 



/ \ -ujcpik 
iuVxpi ■ (j)ik + LoipiV ■ 4>ik] dx 



ib3)k= - 







• 7^ 





-^4>2k 



^IpR 



^i'R 



■K 



■K 







-V • 4,1k 





(^1)*, = 


-L 


A ;• 


K"o*)Ht')-<tO] 


dx 






{A2)t, = 


- L 


Y ° ' 

\^ -UJ(t>lt ^ 


\ -uj4>ij J \ -y ■ (pit J \ -V • (pij J 


(19) 


(^3)*, = 


-I 


( ° ' 






{AA)t, - 


- L 


( ° ' 




dx. 




{A,)t, = 


- L 


( ° ' 


\.7^/^ \ / Y^( Y 




{A,)t, = 


- L 


( ° ' 

A -'^4>2t , 


)-('^0H-v%J-(7)] 


dx. 



dx 
dx 

dx 



-V ■• 



'2 k 



dx 



dx 



- / [ojVxpi ■ (j)2k + wV'/V • 02fc] dx. 



The method for solving for P' and f" can be easily modified to solve for P" and v' . In this case 
the weak equation is 






-uoP" \ -CJS 

-V-v' r'^y -V-T 



dx + I LoT- nP' dS, 
'dr 
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Figure 1. The distribution of the eigenvahies of A and the real parts of the 
eigenvalues oi AI^^A for A^ = 30. 

and all the methods above still apply. In fact, to obtain the new matrix for this formulation, we 
simply change the signs of the blocks A4 and Ag, and the changes in b are mostly reversing signs 
and the roles of the two auxiliary functions ipR and tpj. 

6.1. Other Discretizations. Along with the discretization described in Section |6J we have exper- 
imented with two other implementations in which different basis functions are used to represent the 
variable v. The first of these uses the Raviart-Thomas RT[q] elements described in [5]. We found 
that the resulting finite element matrix is much more poorly scaled, with a condition number ap- 
proximately twice as large as that obtained with the nodal bilinear basis. The second method uses 
the RT^i] elements (also described in [2]). In this case, the higher-order basis functions obviously 
result in a somewhat less-sparse finite element matrix, and the condition number is approximately 
the same as that obtained with the all-bilinear discretization. 



7. Conditioning 



As was mentioned, perhaps the greatest numerical advantage to having a minimization formulation 
for the Helmholtz equation is that the matrix produced by the finite element method is symmetric 
positive definite. This allows for the use of methods such as the conjugate gradient method to solve 
the system. Of course, the use of a preconditioning matrix in the conjugate gradient method can 
speed up the convergence considerably, which is especially important when solving the relatively 
large sparse systems generated by the finite element approach outlined above. 

In our approach, there are three basic types of elements used: bilinear elements, first component 
bilinear vector elements, and second component bilinear vector elements. Each of these types of 
elements interacts with all of the other types, and these interactions are what give rise to the blocks 
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Figure 2. The growth of the number of PCG iterations required to solve a given 
problem with grid size for several error tolerances (outer implementation of PCG 
only). 



in (18 1. Assuming that interactions among similar element types are most important, we choose 
the block Jacobi preconditioner 

/ Ai 
M=i A2 

V ^3 

Among all block diagonal preconditioners of this form, this choice of M minimizes the condition 
number of M^^ AM^^ to within a factor of 3 of its minimum [i]. 

As one of the steps in the preconditioned conjugate gradient method (PCG), [S], a system of the 
form Mr = y must be solved. In order to make solving this problem more efficient, we precondition 
the matrix M and use conjugate gradient to solve this system as well. The preconditioner used 
in this inner implementation of PCG was an incomplete Choleski factorization of M. Figure [l] 
shows the distribution of the eigenvalues of the matrix A before and after preconditioning for 
N — 30. In Figure [2] we see the how the number of PCG iterations grows with N for several error 
tolerances. 

A key component in ensuring that the system Aa = & is well conditioned is for the matrix £ (or 
equivalently TZ and /C) to have a coercivity constant that is as large as possible. For this reason, 
we expect better numerical results when the eigenvalues of C are bounded well away from zero. In 
the case of the Helmholtz equation, the matrix Z is diagonal, say Z = diag(ci, . . . , c^+i), which 
makes it possible to calculate the eigenvalues of C. If D is an invertible matrix, then we may factor 
a block diagonal matrix 

A B 

C D 



as 
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A B\(I B\(A-BD-^C 
C D j ^\ Q D j\ D-^C I 



which imphes that 



det ( ^ n ) =det(L>)det(A-BD-iC). 



C D 



Therefore, 

det(£ - A/) = (-1)''+Mct 



{Z")-^Z' {Z")-^ - XI 

Z" + Z'{Z")-^Z' -\I z'{z")-^ 



= {-lf+^dci{Z'{Z")-^)dct{{Z")-^Z' + [-{Z')-^ + XZ"{Z')-%Z" + Z'{Z")-^Z' - XI]) 
= {~lf+^detiZ'{Z")-^)det{X'^[-Z"iZ')-^] + X[iZ')-^ + Z"{Z')-^Z" + Z'] - {Z'Y^Z"). 
In the case of diagonal Z, this implies that 



where 



-a, ± J a] - h) 

x = — y-^ — - j-i,...,d+i, 



1 fc")^ c" 

a, = - + ^+c; and 6, = 2-f . 



3 3 3 

If Z' = 0, then C is diagonal, and its eigenvalues are those of Z" and {Z")~^. 

The above analysis tells us that the finite element problem will be better conditioned for those 
problems where the coefficients p and k are such that Z is close to li, i.e. p = il and k = —li 
(this would correspond to the limiting case where Uj = bj). In many cases when we are presented 
with a problem where the coercivity constant for £ is small, we can apply an appropriate rotation 
and scaling to the problem in order to get a finite element matrix that is better conditioned. By 
multiplying the problem (l2| through by a complex constant re*^, we effectively replace Z with 
re^^Z, so we should choose r and 9 so that re^^Z is as close as possible to il. However, this may 
not always be possible, for example, when an isotropic p{x) oscillates between values in the upper 
half of the complex plane that are close to 1 and —1. 



8. Numerical Results 

As an example, we demonstrate the error bound on the problem (11]), with parameters p = (—5 + 
5i)/, K = A — Ai and oj = 2. A solution is P{x, y) = e^ii-sy^ j^^ ^j^jg gxample we took 

-0^ = Re(e^""^^) + sin(7ra;) sin(7ry), i'l = Im(e^""^'^) + sin(7rx) • 3sin(7ry)) 

and solved the problem on grids with iV = 3, . . . , 100. Table 1 shows the error in the finite element 
solution for various values of TV. The errors were calculated using the trapezoidal rule with function 
evaluations on a grid with size N = 1500. Figure [3] demonstrates the method on a problem with 
non-constant coefficients, where the dissipation in the material is higher inside a disk centered in 
the unit square. The boundary conditions for the real part are oscillatory, while the boundary 
conditions for the imaginary part are simply an affine function. 
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N 


h 


\\{P-Pn,v-vn)\\v 


30 


0.0345 


6.6162 X 10-^ 


40 


0.0256 


3.6692 X 10"-* 


50 


0.0204 


2.3252 X 10--* 


60 


0.0169 


1.6026 X 10-* 


70 


0.0145 


1.1722 X 10-* 


80 


0.0127 


8.9706 X 10-^ 


90 


0.0112 


7.0686 X 10-^ 


100 


0.0101 


5.7037 X 10-5 



Table 1. The error in the finite element solution for various values of the grid size N. 





Figure 3. The solutions P' and P" , with lo = 10, V'fl — sin(67ra;) cos(37ry), ipj ~ 
3a; + 5y + 2, iV = 30. There is a circular inclusion in the center of the domain 
with p = .01 + .001* and k = .01 — .003i outside the inclusion and p = —5 + 5i and 
K — A — Ai inside the inclusion. 

9. Robin Boundary Conditions 



9.1. Problem Formulation. Another boundary condition that often appears is the Robin prob- 
lem 

-V-p"*VP P = inF, 

K 

P + av ■ n ~ g on dT, 

where a G C In order to deal with this boundary condition, which concerns both real and imaginary 
parts of the variables P and v simultaneously, we start with the minimization functional for the 
natural boundary conditions 

Y{P', v") +2lo ( [P'v' ■ n + P"v" ■ n] dS. 
JdV 

The Euler-Lagrange Equation for the corresponding variational principle is 

B{P\ v", s, T) = -uj [ [sv' ■ n + P"T ■ n] dS. 
Jar 
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where the bihnear form B is defined in (I?]). Notice that we can write the surface integral above 



as 



ar 



V ■ n 
P" 



s 
T -n 



dS. 



The vector on the right contains the primary variables for which we would like to solve, and the 
vector on the left contains the dual variables which we would like to eliminate using the Robin 
boundary condition. In terms of the vectors above, we can express the Robin condition as 



where 



Rearranging, we find that 



Ml 



Ml = 



V ■ n 
P" 



P' 



a' 



= M- 



Mo 



V ■ n 
P" 



9 

g" 



and Mo = 



a' 
a" 1 



- M-^Mi 



P' 

v" ■ n 



so the surface integral term becomes 



— w 



dV 



M2-I 



9 
9" 



M^'^Mi 



P' 



S 
T ■ n 



dS 



dr 



M2-1 



.9 
9" 



P' 

v" ■ n 



dS + Lo I M^^Mi 

IdT 



P' 

v" -n 



P' 



dS. 



The new Euler-Lagrange equation for the Robin boundary condition is therefore 



Jar 



P' 



s 
T -n 



dS 



M^-i 



Since 



we have 



dV 



M^' = 



s 
T -n 



1 

-a" a' 



dS. 



^ a' \ -a" |a|2 

which is positive definite as long as a' > 0. The new bilinear form above is guaranteed to be coercive 
as long as p and k satisfy ( 16 ) and a' < 0. 

To find a numerical solution for the Robin boundary value problem, we discretize using the finite 
element scheme presented in Section |6] Unfortunately, the surface integrals can no longer be 
converted to volume integrals by integration by parts and must be computed as they stand. In this 
case, the finite element matrix is written as the sum of two matrices A — ujB, where A is of the 
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block entries 



form (|18|), and the blocks have entries (|19|), and another matrix B with the same block form and 



(Si)t, 

{B3h 
{B4)t3 
{B5)t3 
{B6)tj 











dS 



dr 



dr 



dr 



dr 



dr 





ht ■ n 


ht ■ n 


ht -n 


ht ■ n 


ht -n 



M^'^Mi 



M^'^Mi 



M^'^Mi 



■ M^^Mi 





(t)ij ■ n 



(t)2j ■ n 





(t>ij ■ n 





dS 



dS 
dS 

dS 



dS. 




The right-hand side vector b is also partitioned as (6i, 62, b^)'^ with entries 

(6l)fe - • ' ' ^« 1 ^r-l I 9' 

{b2)k 
{b3)k 

Assuming that the coercivity requirements ( |16[ ) on p and k are satisfied, and a' < 0, the sys- 
tem 

{A-ujB)a:^b 
may be solved using the same preconditioned conjugate gradient approach as outlined previ- 
ously. 

9.2. Numerical Examples. Here we present some numerical examples obtained by using the finite 
element method to solve problems with Robin boundary conditions. In these examples the Robin 
boundary conditions are imposed on y = and y — 1, while on the other sides of the domain we 
have imposed periodic boundary conditions. On the left in Figure |4] is the solution with a circular 
scatterer with p = 1 + .Ollz outside the scatterer, p — 2 + .Olli inside the scatterer, k = 1 + .Olli 
everywhere, a = —I + .333i and g — 3.33i. On the right, the circular scatterer is replaced by a bar 
angled across the domain, but the other parameters in the problem remain the same. These results 
were calculated using the RTtQi discretization for the v variable described in Section 



6.1 



10. Conclusions 

The variational principles of Milton, Seppecher, Bouchitte, and Willis make it possible to formulate 
the solution of the Helmholtz equation as a minimization, and this is reflected in the fact that 
the stiffness matrix for the finite element method is symmetric positive definite. This allows us to 
use classical iterative methods such as preconditioned conjugate gradient to solve the associated 
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Figure 4. On the left is shown the solution to the Robin problem with a disc of 
high density material centered in the domain. On the right the disk is replaced 
with a bar of the same material angled from the lower left to the upper right of 
the domain. 

system, along with straightforward finite element error estimates. The primary advantage of this 
approach is that it allows the use of efficient iterative methods for the solution of the linear system. 
But there are also disadvantages in that the system has more unknowns, since we must solve for P 
and V simultaneously. 

More research is necessary to determine the circumstances under which this approach may be 
more effective than others currently in use. A particular point of interest is that even though the 
underlying minimization principles are valid for arbitrarily small loss coefficients, the conditioning 
of the associated finite element matrix deteriorates as the system becomes less dissipative. The 
general question of how solution efficiency depends on loss should be analyzed further. 

We have only approached the scalar, two-dimensional Helmholtz equation in this paper, while the 
minimization principles of Milton, Seppecher, Bouchitte, and Willis apply to the full vector Maxwell 
equations, as well as the equations of linear elasticity. The general approach taken here should also 
apply in those cases. We note finally that in many applications for these models, one would like 
to apply nonlocal transmission or radiation boundary conditions in order to accurately handle 
unbounded domains. The problem of adapting these minimization methods to such boundary 
conditions remains open, although presumably the PML approach (see eg. j9|) would apply. 
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